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The Abelian Sandpile Model is a discrete diffusion process denned on graphs (Dhar (19] . Dhar et al. 
£Nj ■ [201 ) which serves as the standard model of self-organized criticality. The transience class of a sandpile 

is defined as the maximum number of particles that can be added without making the system recurrent 
([5]). We demonstrate a large class of sandpile which have polynomially bound transience classes by 
identifying key graph properties that play a role in the rapid diffusion process. These are the volume 
growth parameters, boundary regularity type properties and non-empty interior type constraints. This 
generalizes a previous result by Babai and Gorodezky (SODA 2007, [3]), in which they establish polynomial 
bounds on n x n grid. Indeed the properties we show are based on ideas extracted from their proof as 
well as the continuous analogs in complex analysis. For degree bounded graphs, we demonstrate the first 
constant factor approximation algorithm for the transience class problem, based on harmonic functions. 
This algorithm is used in an essential way in proving the impulse superposition properties, which in turn 
help us derive the polynomial bounds. We conclude with a discussion on the notion of degeneracy and 
dimensions in graphs. 

(N 

^ ■ 1 Introduction 

The abelian sandpile model(ASM) is a type of discrete diffusion process defined on graphs. The model was 
pioneered by Dhar |19| as means of investigating the phenomena known as self- organized criticality, in the 
dynamics of sandpile formation. A close cousin is the celebrated loop-erased random walks model [3D]. If 
one plots the length of a loop-erased random walk, on the same underlying graph, against time, the results 
obtained are qualitatively similar to the one obtained from observing sandpile weights against time. Indeed 
5^ 1 the similarity extends into numerous seemingly different phenomena like stress distribution in earthquakes, 

size distribution in raindrops, etc (see the recent comprehensive survey article by Dhar [18j). 

In the standard sandpile model, "sand particles" are added at vertices of a (multi) graph. A site (vertex) 
is stable as long as the number of particles at the site are less than its degree. Adding any more particles 
makes the site unstable and is accompanied by the unstable site passing a particle along each incident edge to 
its neighboring sites. This relaxation process is called toppling. There is a special site called the sink which 
cannot topple. To ensure that every relaxation process eventually stabilizes, one needs the condition that 
the sink is reachable from every other site. In the course of evolution via particle additions and toppling, the 
sandpile goes through a sequence of configurations. Those which can be revisited in any toppling sequence are 
called recurrent, the remaining ones are termed transient. Typically, one starts with the empty configuration 
and as particles are added, one moves through transient configurations till a recurrent configuration is reached. 
Thereafter the configurations stay recurrent. The steady state behavior of a sandpile is characterized by its 
set of recurrent states. It has been observed by physicists that for most natural phenomena, the time taken 
to reach a recurrent state is small. Hence any acceptable model must reflect this behavior and an important 
parameter of study for these systems is the time taken to reach recurrence. 

The parameter of importance in our discussion will be the number of particles needed to reach recurrence, 
when particles are added adversatively. The strategy of adding particles at uniformly random sites yields 
easily to a simple coupon collector type argument, thereby resulting in polynomial bounds on the expected 
time to recurrence (as previously noted in [5]). However, in the case of adversarial site selections, the problem 
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acquires a distinctly potential theoretic flavor. In this scenario, our goal is to add particles so as to avoid a 
recurrent state for as long as possible. This problem was highlighted by Babai and Toumpakari [5] where they 
define this number of particles as the transience class of the sandpile. This later motivated the insightful work 
by Babai and Gorodezky j3j on grid based sandpiles. In this breakthrough paper, Babai et. al. show that for 
the standard n x n square grid based sandpile, the maximum number of particles one can add before hitting 
a recurrent state is 0(n 30 ), later improved to 0(n 15 ) [2]. Using arguments based on LP-duality, harmonic 
functions, and symmetry, Choure and Vishwanathan [TU] improve this bound to 0(n 7 ). They also establish 
a sandwich theorem for the transience class of any sandpile in terms of the values of harmonic functions, 
showing a tight correspondence between the random walks properties and sandpile diffusions. 

The relevance of grid sandpile comes from the fact that they serve as the standard discrete substitutes of 
planar regions in both statistical physics and traditional engineering research in finite element methods. It 
is therefore natural to ask if other graphs which can arise as regular (or irregular even) tessellations of plane 
(or higher dimensional spaces) also obey polynomial bounds on their transience classes. More generally, one 
might not even want to pose restrictions on the integrality of dimension of the space from which these graphs 
arise. Indeed graphs associated with fractals are quite well analyzed in traditional potential theory. We 
refer the reader to the excellent monograph by Andras Teles (41] for a succinct yet thorough introduction to 
literature in this area. The techniques of Babai and Gorodezky [3] extend easily to the honeycomb lattices 
and triangular lattices. They can also easily be generalized to higher dimension lattices. The problem arises 
when one wants to deal with graphs without such highly symmetric structure. Our study is aimed at trying 
to understand the properties which characterize polynomial transience. 

Our contribution: Using LP duality and the independent set properties of nodes with zero heights, 
we demonstrate an algorithm which approximates the transience class of degree bounded sandpile up to 
constant factors. The algorithm works by computing harmonic functions over graphs and is the first such 
result. This concludes the status of our previous bound on transience classes |10| and significantly tightens 
the connection between sandpile and random walks on the underlying graph. We show that if the graph has 
uniform polynomial volume growth, satisfies mean value property, and has high local conductance of particle 
percolation, then it satisfies a fundamental superposition property, which is the discrete analogue of the 
celebrated superposition principal. This property allows one to place particles at a suitable set of multiple 
nodes and observe the same potential response at a particular site, as one would have by placing all these 
particles at one site. Using this we prove one of the main results in this paper, the single step Epicenter 
Propagation Lemma. This is the general form of the traveling diamond lemma proved in Babai and Gorodezky 
[3]. This lemma demonstrates a very general gluing property of potential functions. No continuous analogue is 
known for harmonic functions over continuous regions, where to bound the growth rate of harmonic functions 
over any path, the best one can do is simply multiply the bounds over each of the balls forming a cover of the 
path. Using the single step version, we derive the general Epicenter Propagation Lemma which effectively 
demonstrates polynomial bounds on the transience class of large class of graphs. This result forms the graph 
theoretic analogue of the classical Harnack's inequality and enables one to prove general bounds on harmonic 
functions without resorting to limiting arguments linking random walks on graphs to Brownian motions on 
Euclidean spaces. In this sense, the applicability of our results is more general as there exist lot of interesting 
graphs, for example the ones arising from fractals like the Sierpinski gasket, which cannot be likened to a 
Euclidean space of any integer dimension. We conclude with a discussion on two very fundamental questions 
related to the dimension theory of graphs. The first is to determine a dimension like quality of graph which 
is consistent with the more usual notion of dimension for Euclidean spaces. The second concerns determining 
if a subgraph is degenerate with respect to the ambient super-graph. In the last section, we discuss some 
open questions in this area. 

1.1 Related Work 

Random walks and Sandpile: The connections between sandpile and random walks have been discussed 
under various contexts. Dhar, for example, summarizes the connection between the Loop Erased Random 
Walks and sandpile in [18]. We refer the reader to the excellent text by Lawler and Limic [30] for an in- 
troduction to basic properties of LERW. The connection with simple random walks is much more intuitive 
and reveals a beautiful analogy between properties of sandpile and electric networks. The recent results by 
Choure and Vishwanathan ITU] discuss a reduction of the transience class computation problem to that of 
estimating potentials in a related circuit. These connections are not surprising as the combinatorial Laplacian 
of the graph of underlying sandpile acts as kind of generator for both the models. In the present paper we 
strengthen this correspondence by presenting the first constant factor approximation algorithm for estimating 
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transience class sizes in terms of potentials, for bounded degree graphs. For an introduction to these potential 
functions, see Bollobas [9]. For a discussion on estimating these potentials on graphs coming from geometric 
settings, the monograph by Teles [41] is highly recommended. 

Diffusions and Potential Theory on Graphs: Delmotte [UJ shows how Harnack's inequality leads to 
non-trivial results on heat diffusions. They talk about bounding the growth of harmonic functions using 
Gaussian estimates of the heat kernel, all of which follows from the assumption that the graph under con- 
sideration follows a parabolic Harnack's inequality. Chung and Yau |12| derive the Harnack inequality for 
certain degree bound graphs. The bounds are then used to bound the Neumann eigenvalues of these graphs. 
In [TS], they derive Harnack's inequality of Abelian homogenous graphs. Further more, Chung and Yau [T4] 
derive lower bounds on log-Sobolev constants by establishing log-Harnack inequalities on graphs. Bounding 
the log-Sobolev constants is important as it helps in establishing better convergence bounds on random walks 
on graphs, see Diaconis and Saloff-Coste [5T]. Chung and Yau [TJJj derive bounds on the eigenvalues of the 
Laplacians using the Sobolev inequalities and heat kernel estimates. 

Electric Networks: The classical theory of random walks ([32], [22], [31]) has some very powerful and intu- 
itive results which have recently found widespread application in theoretical computer science. Christiano, 
Kelner, Madry and Spielman |11| have recently announced the fastest known algorithm for computing approx- 
imate maximum s — t flows in capacitated undirected graphs. Their algorithm constructs the approximate 
flows by using the electric current flows on the same network with s and t as poles in an essential way. Earlier 
Kelner and Madry [26] used arguments based on random walks to formulate the fastest known algorithm for 
generating spanning trees from uniform distribution. Spielman and Srivastava |38) construct good sparsifiers 
of weighted graphs via an efficient algorithm for computing approximate effective resistance between any two 
vertices, a result which is quite insightful on its own. Indeed a deeper understanding of harmonic functions is 
as much of interest to a computer scientist as to a potential theorist and as such the benefit of this confluence 
of different types of researchers has been mutual. A significant example is the, now seminal, work by Arora, 
Rao and Vazirani [1] on embeddings of negative type which give an 0(^/\ogn) approximation algorithm for 
computing graph conductance. 

Other Results in Sandpiles: As already mentioned, research problems in abelian sandpile model span 
across numerous areas. We make a passing mention to some of this work. Recent advances in complexity 
theoretic flavor include proof of the one-dimensional sandpiles prediction problem in LOGDCFL by Peter 
Bro Milterson [34]. Schulz [36] mentions a related NP-complete problem. The group structure of the space 
of recurrent configurations, first introduced by Dhar, Ruelle, Sen and Verma in [20] is also considered a 
fertile area of analysis. Cori and Rossin |16| show that sandpile groups of dual planar graphs are isomorphic. 
Toumpakari [43] discusses some interesting properties of sandpile groups of regular trees. Questions related 
to group rank are studied in particular, the paper is concluded with an interesting conjecture on the rank of 
all Sylow subgroups of the sandpile group. Specific families of graphs like square cycles C^, K$ x C n , 3 x n 
twisted bracelets, etc have been analyzed. We refer the reader to [25], [37], [H]- 

2 Basic Properties of The Abelian Sandpile Model 

Our notation and terminology follows Babai and Gorodezky [3J. 

Definition 1. A graph G is an ordered pair (V(G), E(G)) where V(G) is called the set of vertices and E(G) 
is a set of 2— subsets of V, possibly with repeated elements, the set of edges. 

This is referred to as a multi-graph in literature but we will use graph for brevity. The degree of a vertex 
v G V is defined as the number of edges in E which contain v. Two vertices v and u are called adjacent (or 
neighboring) if (it, v) G E. A path between two vertices u and v is an ordered sequence of edges ei, 62, . . . , &k 
such that u G e±, v G and for all values of i, e, l~l e i+1 ^ <f>. The graph G is connected if there exists a path 
between any pair of vertices. 

For an Abelian Sandpile Model, we take a connected graph G with a special vertex called the sink, denoted 
s G V . Non-sink vertices in G are called ordinary vertices and this set will be denoted by V = V — {s}. 

Definition 2. The configuration of a sandpile G is a map c : V Z + . It will be represented as a vector. 
The weight of c is \c\ = J2vev c ( v )- 
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The configuration c tells us the number of sand particles that each of the ordinary sites currently contain. 
The empty configuration is the zero vector. The capacity of a site is the maximum number of particles that 
it can hold and is one less then the degree of the node. 

Definition 3. An ordinary node v is said to be unstable in a configuration c if c(v) > degree (?;). If all the 
sites in a configuration stable, the configuration is stable, else it is referred to as unstable. 

When a site is unstable it is said to topple, i.e. pass on some of its particles to its neighbors. When a site 
v topples once, it loses degree(u) particles and each neighbor of v acquires a particle for every edge common 
with v. The sink node, by definition, never topples. We start with the empty configuration and keep adding 
particles one by one on sites of our choice and topple when necessary. 

The ASM evolves in time through two modes, particle addition at sites and relaxation of unstable sites 
via topplings. A toppling sequence is an ordered set of configurations where every configuration can be 
obtained from the previous one by toppling some site unstable in it. The case of many sites becoming 
unstable simultaneously also poses no complication as the order in which they are subsequently relaxed does 
not effect the final stable configuration that is obtained at the end of toppling sequence, hence the prefix 
abelian. Elementary proofs of such confluence properties can be found in the pioneering paper on ASMs by 
Dhar [IH]. 

Notation: We write c\ > C2 if Vu, Ci(t>) > ci and c± h C2 if there is a toppling sequence which takes C\ 
to C2- Lastly we write, c\ — > ci if 3c3 > c\ such that C3 h C2. We say that a configuration C2 is reachable 
from ci if c\ —¥ C2 and unreachable otherwise. In words, one can add particles to certain sites in ci so 
that there exists a toppling sequence leading to C2. Note that being reachable is a transitive relation, i.e. 
ci -> c 2 , c 2 -> c 3 => ci -> c 3 . 

Theorem 2.1. (\Wj,JBj) Given any configuration c, there exists a unique stable configuration o~(c) such that 
c h cr(c), independent of the toppling sequence chosen. 

Property 2.1. If c h cr(c), then kc h ka(c) 

Property 1 2. II will later be extensively used for proving gluing properties of potential functions. Associated 
with every toppling sequence is the count of the number of times each site has toppled, the vector of toppling 
potentials, also referred to as the score vector in [3J. These toppling potentials are very closely related to 
the electric potentials that develop at various nodes when power source-sink are appropriately applied, a 
connection which we will discuss in detail in the coming sections. 

Definition 4. Assuming c\ h C2, the toppling potential function z Cl ' C2 : Vq —> Z + is defined as z Cl,C2 (v) : the 
number of times v toppled in a toppling sequence from c\ to C2. We denote z c,a ^ by z c . 

This function is well defined as the number of times a particular site topples is independent of the toppling 
sequence chosen, already noted in [3j . A simple proof follows from typical linear algebraic arguments and the 
fact that the principal minor of a connected graph's combinatorial Laplacian is of full rank. 

A configuration is called recurrent if it is reachable from any configuration. As already mentioned, 
we say that a configuration Cj is reachable from a configuration Cj if by adding some particles to Cj and 
subsequently relaxing it, we can obtain Cj. A configuration is transient if it is not recurrent. The set of 
recurrent configurations is therefore, closed under being reachable. 

Definition 5. A configuration c is recurrent iff Vc' we have c' — > c. 

Since recurrence persists under particle addition, we have the following property. 

Property 2.2. If configuration c\ < C2, then recurrence of c\ implies that of ci. 

Denote the configuration in which every node v has degree(u) particles by c max . Clearly, given any stable 
configuration c, one can reach c max simply by adding the required number of particles at each site. 

We analyze the process of adding one grain at a time to the sandpile and study its evolution. As in 
the standard theory of Markov chains, recurrence characterizes the long term (steady state) behavior of 
sandpile. Our investigation is concerned with the maximum number of particles that can be added while 
staying transient. Following Babai and Toumpakari [5], for a sandpile S we define the notion of transience 
class as follows. 

Definition 6. The transience class of S denoted by tcl(S'), is defined as the maximum number of particles 
that can be added to S before reaching a recurrent configuration. 
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Remark: In |10| , the authors talk about an alternate equivalent characterization of transience class which 
defines the transience class alternatively as the maximum number of particles that can be added before all the 
nodes have toppled at least once. Even though we will not be using this definition explicitly, the notion is 
inherent in the way we bound transience classes. We will be computing the maximum number of particles 
that can be added at any point before at least one particle reaches every node. For sandpile with bounded 
node degrees, using the Property (|2.ip . this translates to the alternate characterization of transience classes 
we just mentioned (albeit with a constant multiplicative loss factor, the degree). 

3 Basic Potential Theory 

For a proper introduction to harmonic functions on graphs, we refer the reader to the beautiful paper by 
Benjamini and Lovasz [6|. See Teles |41| for a thorough view. We start with some important definitions and 
fundamental properties. Given a connected graph G and a function tt : V(G) — > R, we say that ir is harmonic 
over the vertex set V/j if, 



The vertices in V — Vh, adjacent to any vertex in Vh, are called the "poles" of ir. The set Vh is also called 
the interior of 7r and the set of poles referred to as the boundary. Being harmonic over Vh means that the 
value of 7r at any vertex in Vh is the average of its value in the immediate neighborhood. In case of multi 
graphs, we take the appropriate weighted means, where the weights are the number of common edges. This 
leads us to the first basic property, 

Property 3.1. Any non-constant harmonic function can assume its extreme values only at the set of poles. 

It follows that every non-constant harmonic function has at least two poles, its maxima and minima. 
Such functions are completely determined by their values on these vertices. Formally speaking, 

Property 3.2. Uniqueness: If two functions harmonic on Vh agree on the boundary, they agree everywhere 
in the interior. 

Property 13.21 is important as it allows one considerable freedom in constructing harmonic completions of 
functions defined over the boundary set. This problem is the discrete analogue of the classical boundary 
value problems in complex analysis. If one fixes the interior set Vh and allows arbitrary boundary values, the 
complete set of harmonic functions is obtained. This set is closed under linear combinations and contains all 
the constant functions. 

Property 3.3. The set of functions harmonic over any set Vh form a vector space. 

Two important scenarios in which these functions arise naturally are electric networks and random walks. 

Electric Networks: Consider a resistive electric network (i.e. a circuit made up entirely of resistors). Let 
s TTt(v) be the potential that appears at node v when unit potential is applied across t and s. Using the 
equation of charge conservation (Kirchoff's node law), one can show that these potentials are harmonic on 
all nodes except s and t. 

Random Walks on Graphs- dipole version: Consider a graph G and two special vertices s and t. The 
potential associated with v, with s and t as poles, s 7r t (w) is defined as the probability of reaching t before s 
starting from v. One can check that the function 7r so defined is indeed harmonic on the set V — {s, t}, with 
the maximum value of 1 at the node t and the minimum value at s. The generalization to the multi-pole 
situation is simple, but it is interesting for a different reason. 

Random Walks on Graphs- multipole version: Again we consider graph G, but this time we have a set of 
interior nodes Vh and the corresponding set of poles, S = {sj}. Denote by iTj(v), the probability of a random 
walk, starting at v, hitting Sj before any other vertex in S. Assign any desired set of values 7r(sj) to each of 
these poles. The value ir(v) is defined as the sum ir(sj)TVj(v). As before, it is a simple exercise to check 
the function n is harmonic over Vh . The formula basically suggests taking an expectation over the boundary 
values, using the measure induced by the random walk. 

The main implication here is that one can intuitively think of the electric network theory as an analysis 
of random walks of electrons on the underlying graphs. Consequently, results from network theory can be 
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used to prove interesting facts in other related areas. As an example, consider the problem of constructing 
the harmonic completion of a function with given boundary values. All one needs to do is to take the 
corresponding circuit and apply potentials equal to the boundary values on the boundary points. The 
potentials that will appear on other nodes can be computed using basic linear algebra (the only non-trivial 
step involves inverting the combinatorial Laplacian of G) thus allowing construction of harmonic completions 
efficiently. 

The reciprocity theorem can be restated in terms of just potential sources and potential measurements 
using the notion of effective resistances between pairs of nodes. The effective resistance between a pair of 
nodes u and v, R e ff(u,v) is defined as the potential difference which develops between u and v if a unit 
current source is applied across u and v. In any resistive network, the following reciprocity property holds. 
See [TU] for an elementary proof. 

Lemma 3.1. Potential Reciprocity Lemma : If taking s and t as poles with 7r(s) = and ir(t) = 1 induces 
a potential of s TTt(v) at node v and interchanging the roles of v and t induces s Tr v (t) at t then, 

Reff{s,t) s TT t {v) — R e ff(s, v) s TT v (t) (1) 

In particular, when the effective resistances across s and t are the same as s and v, we have s irt(v) =s Tr v (t). 
In the following discussion, we will omit the left subscript (s) from s ir t whenever it is clear from context. We 
say that a walk P is an instance of s ir t if it starts at some vertex v, avoids s and ends at t. The following 
lemma may already be known to experts. A proof appears in [10] . 

Lemma 3.2. A triangle inequality for potentials 

Ki{j).Kj{k) < 7T<(fc) (2) 

4 Basic Properties 

Consider any (possibly infinite) graph G(V,E). Distances in this graph will henceforth correspond to the 
shortest path metric. For any vertex v, B(v, r) denotes the ball of radius r around v, i.e. the set of all vertices 
in G which are at a (shortest path) distance of at most r from v. For any set of vertices U C V, define two 
notions of boundary 

- vertex-boundary, Syll is the set of vertices, in U, having neighbors in the set V — U . 

- edge-boundary, SeU is the set of edges connecting vertices in U to those in V — U. 

For example if the set U is the ball B(v,r), the set 6vB(v,r) is the set of vertices which are exactly at 
distance r from v and the set SeB(v, r) is the set of edges between the vertices in 6vB(v,r) and those in 
SyB(v,r + 1). Consider any connected set of vertices Vh and the induced subgraph Gh = (Vh,Eh). The 
sandpile corresponding to this sub-graph is obtained by adding the edge set SeVh to Eh and identifying all 
the vertices in V — Vh which have neighborhood in Vh with the sink node s. We denote this sandpile by 
S=(V h U{s},E h U5 E Vh). 

Example of grid sandpile: To make the above definition more transparent, consider the example of grid 
sandpile. Consider the infinite grid graph with its canonical embedding in plane. For every pair of integers 
(i,j) their exists a vertex (with this pair as its label). Each vertex is adjacent with (i + (i,j + 1), 
(i — 1, j) and (i,j — 1), which are just the four lattice points flanking Let V m ^ n be the set of vertices 

with labels ■ < i < to,0 < j < n}. Evidently this set is connected. The boundary vertex set is the 

set of vertices lying on the horizontal line segments {(i, 0) : < i < m} and {(i, 0) : < i < m} and on the 
vertical ones {(0, j) : < j < n} and {(m,j) : < j < n}. The edge-boundary set is the set of edges between 
vertices in the boundary set and the vertices not in the (m + 1) x (n + 1) block (Vh). Figure Q] exhibits 
the graph, the aforementioned boundary sets and the resulting sandpile. In the left figure, the white nodes 
belong to the set V — V m , n . The grey nodes form the vertex boundary of V m . n , and along with the black 
nodes constitute the set V m . n . The heavy edges constitute the edge boundary. The figure on the right is the 
sandpile S m , n - The heavy edges are connections to the special sink node. 

We will denote both the sandpile and the underlying graph by the same symbol S. To keep notation 
and definitions clean, we fix the following convention. A ball B(v,r) inside G is defined already. However, 
when we say a ball B(v,r) in some vertex subset Vh, or the induced sandpile S, we assume it is small 
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Figure 1: A grid graph, the subset set V m ^ n and the corresponding sandpile 



enough to not contain the sink node. Also, whenever a sandpile is mentioned without explicit description 
of an ambient graph, we assume that a super-graph indeed exists and the sandpile has been obtained in the 
manner described above. We now define some basic properties which we will need later in discussion on 
bounding the transience class. 

Definition 7. Polynomial Volume Growth Property (V a ) : A graph G has polynomial volume growth 
property (V^,) with growth exponent a if there exist constants 8 and A such that for any node v and ball 
B(v 1 r) inside it, 



S.r a < Vol(B(v,r)) < A.r a 
Here, Vo\(B(v, r)) is the number of edges inside the ball B(v, r) 



(3) 



This property plays two crucial roles. First, it ensures measure isotropy in a graph. Which means 
informally that no matter where we are in the graph, if we start a random walk, it can go in any direction 
more or less uniformly. So there is no drift in any particular direction. Secondly, there is a bound on the rate 
at which the measure grows. This is crucial if one wants to make any arguments along the lines of diffusions 
in Euclidean spaces. An easy consequence is the degree bound on graph. All vertices have degrees bound 
between S and A. 

Definition 8. Degree boundedness (A) : A graph satisfies (A) if every node's degree is bounded from above 
by A. A sandpile satisfies (A) if every normal (non-sink) node has degree bounded by A. 

If adding x particles at some node u causes some node v to receive a particle at any point in time, we say 
that v is flooded. Similarly, we say that the set S C V got flooded if it consists of sites all of which received 
at least one particle. There are two symmetry properties that will be important to us. 

Definition 9. High Local-Conductance Property (hLC) : A sandpile S has (hLC) if there exists a constant 
C s such that for any site v and any ball B(v,r) in S, placing at most C s .Vo\(B(v, r)) particles at v floods 
the ball B(v,r). 

This property implies that there is not too much dissipation of particles through the boundary of the 
sandpile, as long as we restrict ourselves to ensuring flooding balls that lie completely inside the set of ordinary 
vertices. For the second property, consider any such ball. If one starts a random walk at the center, the 
probability of hitting any vertex on the boundary should not be much different from any other. This is quite 
like the measure isotropy type property we have previously mentioned. For our proofs, we will be needing 
the weaker form of this property. 
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Definition 10. Mean Value Property (MV) : A sandpile graph 5" = (V h U {s}, E h U 8 E V h ) has (MV) if there 
exists a constant C such that for any function tt : V — > [0, 1] defined over V which is harmonic over Vh, the 
following inequality holds for every ball B(v,r) inside Vh, 

ir(u)>C.n(v).Vo\(B(v,r)) (4) 

The continuous analog of mean value properties occupy central position in classical analysis. See Section 
ijSJ) for further discussion on the importance of this assumption. We now define an impulse superposition 
property. 

Definition 11. Local Superposition Property (LS) : Given a ball B(v,R) C S, if placing H particles at v 
topples a site w in B(v,R) and placing h particles at every site in the smaller ball B(v,r) also topples w, 
with h minimum, then there exists a constant C such that 

h * a wk^) (5) 

Remark: This property essentially means that if placing H particles at a site v causes a toppling at site w, 
then the same effect can be produced at w by adding a constant factor times H particles distributed evenly 
among the sites in the ball B(v,r). Such impulse distribution properties are easy to prove in the setting of 
classical complex analysis. Our more ambitious goal of proving a very general superposition property where 
the impulse is distributed not over a ball but a general set of nodes which are distance wise (in potential 
theoretic sense) well distributed with respect to the point of observation, did not work out. It seems that the 
constraint of taking nodes inside a ball is essential. 

Theorem 4.1. MV + A ^>LS. 

One would expect this property to hold if one considers the sandpile as a discrete model of random walks, 
as shown in |10] . However it is surprising in the light that this property is a hallmark of linear systems, 
a family of which random walks based models are but illustrious members, while sandpiles are models of 
non-linear dynamical systems [18] . We defer the proof of this lemma till the constant factor approximation 
algorithm for transience classes is discussed. 

Definition 12. Overlapping Potentials Property (OP) : Given a sandpile S, let h be the smallest number of 
particles which when placed at every site in a ball B(V, r) C S makes every site in the normal ball B{v, R) C S 
topple at least once. If h is bounded by a polynomial in R/r, the sandpile S satisfies the overlapping potentials 
property. 

Lemma 4.1. hLC + MV + V a ^OP. 

Proof: Assume the sandpile satisfies hLC, MV, and V a . Using theorem 14. 11 it also satisfies LS. Consider 
any vertex v with the balls B(v,r) and B(v,R) in S. Let x be the minimum number of particles one needs 
to place at v to topple every site in B(v, R) at least once. Consider any site p in SB(v, R). 

Using (LS), if x particles placed at v induce a toppling at p, then (f> s .x/Vol(B(v, r)) particles at each of 
the sites in B(v, r) necessarily induce a toppling at p. Since S satisfies (hLC), we have 

x w C s .R a 

Similarly, the property (V^) implies 

Vol{B(v,r)) w r a 

So the number of particles needed to be placed at each node of B(v, r) to topple every site in B(v, R) at 
least once is, 

w 

Therefore, h is indeed a polynomial in the ratio R/r. □ 
Notation: With each site v in S we associate the number r){v) = dist(i>, Sv(S)) i.e. the distance between site 
v and the vertex boundary of S. It is also the radius of the largest ball centered at v which is inside S. 



Lemma 4.2. (Epicenter Propagation : Single Step) Let v be a node in S which satisfies (OP), with 
B(v, f(v)) as the largest ball around it. Let u be any site in SvB(v, \j)(v)/2)\ (i.e. at distance \j](v)/2\ from 
v). Then, if a configuration c floods the ball B(y,r)(v)), there is some constant K such that the configuration 
K.c floods B(u, rj(u)). 

Proof: To keep the notation clean, we drop the floors. The argument without this simplification needs 
no new ideas and can be essentially reconstructed from the following proof. 

Consider any shortest path connecting v to u. The value of rj can vary by at most one with every 
step. Therefore, rj(v)/2 < 77(14) < 3n(v)/2. Given that the configuration c floods the ball B(v,r/(v)). Since 
B(u,n(v)/2) C B(v,r/(v)), c also floods B(u,rj(v)/2). Using the overlapping potentials property (OP), we 
know that the number of particles needed at every site in the ball B(u,n(v)/2), to flood the normal ball 
B{u,rj{u)), is bounded by a polynomial in 2r/(u)/r](v). Since this ratio is bounded by 3, the polynomial itself 
is bounded by a constant, say K. Using Property 12. 11 K.c places K particles on each site of B(u,rj(v)/2), 
which in turn necessarily floods B(u,i](u)). □ 

Note that there are two distinct case here, n(u) < r/(v) and n(u) > r/(v). Even though one argument 
suffices to obtain bounds for both the cases, they are qualitatively different. See figure ©. The case in which 
the distance from sink node is increasing along the path is termed the expansion step and the other case is 
called the contraction step. There is a third variant also, a drift step. The need for such distinction will 
become clear when we prove the general multi-step version of this lemma in the next section. 



f(v)/2 



f(u) 











i > 




• f(v)/2 f(u) , 









(a) (b) 

Figure 2: Epicenter Propagation : Single step (a) Expansion (b) Contraction 

Remark: The above lemma allows one to shift the focal point of diffusion by some distance at the cost of an 
additional multiplicative constant. This generalizes the traveling diamond lemma of Babai and Gorodezky [3]. 
This lemma forms the backbone of the main result. While proving the polynomial bounds on the transience 
class, we will be using this lemma to keep shifting the epicenter from some starting node to a target node. 
We will bound the total number of applications of this lemma by 0(logn) thereby bounding the total particle 
requirement in poly(n). 



5 Bounding the transience class 

In the proof of the Lemma (|4.2p we use the fact that along any shortest path joining a vertex v to a vertex 
in B(v,r) (lying inside the sandpile), the distance from sink cannot increase or decrease too fast. For our 
purpose, this precondition of u belonging to B(v, r) is too restrictive. Instead, we simply define the following 
property for a path in graph and using it, impose a more general constraint on the graph structure. 

Definition 13. The function 77 : V — > M varies linearly over path P (from v\ to V2) if for any site v at a 
distance distp (ui , v) along P from v\, we have the following two sided bounds on n(v). 

ai + v( v i) + 6-distp(t>i, v) < T](v) < a u + r)(v\) + b.dist P (v 1: v) 
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Here a u and ai are constants with absolute values bounded by 1. 

Note: The value of b lies in the closed interval [—1, 1] as the distance from sink cannot increase (or decrease) 
by more then one, when one step is taken along the path. 

Call a path P from v\ to v-i central if i](vi) values vary linearly over P, or if its length is bounded by 
0(log(\S\). We call it k— central if it is a juxtaposition of at most k central paths, where k is a constant 
independent of graph size. 

Definition 14. Non-empty Interior Property (NI) : A sandpile S has (NI) if for every pair of vertices v and 
w, there exists a path between them which is k— central. 

Definition 15. Polynomial Transience Class Property (pTcl) : A sandpile S satisfies (pTcl) if by adding at 
most polynomial (in sandpile volume) particles at any site, one can induce a toppling at any other site. 

The following lemma proves the polynomial bound on transience class of certain sandpile families. It is 
a generalization of the expansion and contraction phases while flooding a grid, as mentioned in Babai and 
Gorodezky [5]. 

Lemma 5.1. (Epicenter Propagation lemma) OP + NI^> pTcl, i.e. existence of k- central paths between 
every pair of vertices in S along with the overlapping potentials property imply (pTcl). 

Proof : First consider the case k = 1 . The general case will follow immediately by taking the k th power 
of particles required in the special case. Consider the sandpile S and a pair of vertices p and q with a central 
path P from p to q. Denote the vertices on the path P as vq = p,vi, . . . ,v n = q. The value of / changes 
linearly from Tj{p) = a to 77(g) « a + b.n. There are three essentially different cases which can arise, b > 
and b < 0, when rj varies linearly, and path length is logarithmic in graph size, when r) varies sub-linearly. 
We consider them individually. 

1. Case (b > 0) Expansion Phase : Start by adding some constant number of particles so that the normal 
ball B(v 2 ,2) is flooded. Let this constant be K s . We will now repeat the following process iteratively 
till q receives a particle. Assume that we are currently at vertex Wj and the normal ball B(vi,r/(vi)) 
is flooded. Consider the vertex Vj with j = i + n(vi)/2 and the normal balls B s = B(Vj,r)(vi)/2) and 
B t = B(vj,r](vj)). B s was flooded in the last iteration. We use this to flood the the concentric ball 
B t . Using Lemma [4.2i if configuration c flooded B(vi,r](vi)), then K.c floods B(vj,n(vj)). In every 
iteration, the radius of the target ball B t increases by at least a factor of n(vj)/r)(vi) ~ b. If path length 
is n, the total number iterations is at most log(n). Also, each step contributes a bounded multiplicative 
factor to the total count of particles needed. This implies that the total number of particles needed is 
of order K s .K l °o( n ) ~ poly(n). See figure ©. 

2. Case (b < 0) Contraction Phase : This case is analogous to the Expansion phase, except that we add 
0(|5|) particles in the first step to flood B(p, r/(p)) (follows form the (hLC) property). There is constant 
factor contraction in each iteration. The implication of this remains the same though, that the total 
number of iterations are logarithmic in n. Again, using the Lemma I4.2[ each iteration contributes at 
most a constant multiplicative factor to the total particle count. Consequently the total number of 
particles required to achieve complete contraction is poly(n). See figure 

3. Case (|P| = 0(log(n)) Lateral Drift Phase: Let the path P be such that r\ varies sub-linearly over it. 
Let the minimum value of r\ over P be r\ m in- Consider the ball B{p 1 r] m in)- We flood this using at most 
Ko = 0(|5|) particles (follows from the (hLC) property). If configuration c floods B(p, r] m i n ), then K.c 
floods B(^ mjti / 2 , T) m i n ), using Lemma 14.21 At a multiplicative cost of K, we cover a distance of r) min /2 
on the path P. Since \P\ is bounded by 0(log(n)), we take 0(log(n)) steps in this fashion. The total 
number of particle required to finally flood the site q is 0(n.K 0< - l ° 9< ' n ^) which is poly(n). See figure 
©• 

After considering these cases of central paths, the argument of fc-central paths follows intuitively. We 
split the path in the central components. Each of these components contributes a polynomial factor in the 
total particles count. Since k is a constant for the family of sandpile under consideration, we obtain that the 
total particle count is also polynomial. 

□ 
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q 



Figure 3: Epicenter Propagation : Expansion Phase - the radii of balls is expanding and the potential focus 
is moving towards the center of graph interior 




Figure 4: Epicenter Propagation : Contraction Phase - the radii of balls is going down and the potential 
focus is moving from the center of graph interior towards the boundary region 

6 A constant factor approximation for transience classes of degree 
bounded graphs 

In this section we will derive a constant factor approximation algorithm for single site particle addition 
transience classes. We later generalize it to multi site particle addition strategies. Using this, we will derive 
conditions required for a graph to satisfy the local superposition property. 

Consider a sandpile S = (Vh U {s},Eh U SeVh) with the sink node, s. If x v is the maximum number of 
particles that can be added at node v without toppling the node w, then the following linear program bounds 
the value of x from above. 

max x v 

< Z W) — degree(w).z(u) + x v < degree(v) — 1 

\fu 7^ v : < z {u') — degree(u).z(w) < degree(u) — 1 
z(w) < 0,z > 0,x v > 

Note that, in the program above and the ones to follow, the variables v,w, . . . all denote nodes from Vh- 
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Graph Boundary 




Figure 5: Epicenter Propagation : Lateral Drift Phase - the radii of balls under consideration stays same but 
the potential focus drifts inside the central region of the graph from one focal point to another 

Similarly the summations indexed by nodes run over the vertices from Vh only. The node s is special, and 
whenever it figures, we will mention it explicitly. With the fixed sink node, s, we define tt w (v) as the potential 
at node v when a unit potential is applied (directly) at node w. In terms of these potential functions, the 
following bound holds. We give the proof here for completeness. 

Lemma 6.1. (J 101/) The optimum value of above LP is bounded from above by the following value, 

— "TT ^2{degree(u) - l).ir w (u) (7) 
n w {v) u 

Proof: We consider the following relaxed version of the given LP. 

max x v 

z(v') — degree(v).z(v) + x v < degree(w) — 1 

Vit 7^ v : z(u r ) — degree(u).z(u) < degree(u) — 1 
z(w) < 0,z > 0,x v > 

From the weak duality for LPs, it follows that to obtain an upper bound of a on the optimum value of 
the above system, it suffices to find a feasible solution of the dual LP of value a. The dual is the following : 

min y^(degree(tt) — l).Y(u) 

u 

Y(u')+Y' -degree(w).Y{w) > 
Vu^w: V Y(u') - degree(u).Y~(ti) > 
Y» > 1,Y > 0,Y' > 

Consider the following set of equations 

Y Y(u')+Y' -degree(w).Y(w) = (8) 

It,' r^W 

Vu^w: ^ F(u') -degree(u).r(u) = (9) 

u' ~u 

Y(v) = 1 

A non-negative set of values satisfying the above set is feasible for the dual LP. We find these by considering 
the resistive circuit S, obtained by replacing each edge in S by a unit resistance. We assign ground potential 
to the sink, and inject current at node w such that it gets unit potential. The potential that develops on 
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any node u is s n w (u). The potential value at node v, n w (v), can be used to scale the input current at w 
thereby scaling all the potentials as well, such that potential at node v becomes unit. It follows that the 
values Y{u) = tv w {u)/tt w (v) and Y' equaling the value of the current injected form a feasible solution of the 
dual LP. The objective value at this point is, 

— TT y](degree(w) - 1).^^) 
ir w (v) ^ 

□ 

We will now argue that in the case of degree bound graphs, using some combinatorial properties of 
independent sets in sandpile graphs along with arguments similar to the ones outlined above for upper 
bound, one can derive tight (up to constant factors) lower bounds for the transience classes. 

In the sandpile S, let x' v be the minimum number of particles that can need to be added at node v such 
that the node w topples at least once. Clearly x' v = x v + 1. As in [TU], the following LP relaxation forms a 
lower bound on the value of xL. 



< z ( v ') ~ degree(w).z(u) + x v < degree(u) — 1 
Vu v : < z ( u ') — degree(w).z(w) < degree(u) — 1 
z(w) > l,z> 0,x' v > 

Consider the scenario in which we add x' v particles at v and allow the configuration to settle down 
to stability. Furthermore, let the number of particles that appear at any node u in the resulting stable 
configuration be h u . The following lemma gives a lower bound on the particle count x' v in terms of these 
height functions h u . 

Lemma 6.2. The value of the particle count x' v is bounded from below by the following value, 

-J—^/^).^) (10) 

Proof: Clearly < h u < degree(w). In the linear program stated above, replacing each pair of constraints 
of the type < Yt, v '~v z ( v ') ~ d( v )- z ( v ) + x v < d(v) - 1 by J2 V >~ V Z W) ~ d(v).z(v) + x v = h v , maintains 
the property that the optimum value is at least a lower bound to the exact solution to the particle count x' v . 
The altered LP is the following, 



z{v') — degree(w).z(w) + x' v = h(v) 
Vu 7^ v : z(u') — degree(u).z(u) = h(u) 

u' ~u 

z(w) >l,z>0,x' v >0 

To bound the optimum value, we will consider the dual of this minimization program and find a suitable 
feasible point. The value of the cost function at that point will be used as bound. The dual is the following 
maximization program. 

max7' + ^%).y(ii) 

U 

Y{u) + Y' - degree(u;).r(u;) < 

Muj^w: ^ Y(v!) - degree(u).y(w) < 
Y(v) < 1,Y > 0,Y' > 
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As in the previous lemma, we consider the following set of equations, whose feasible region is inside the 
one corresponding to the dual we mentioned above. 

Y(u') + Y' - d{w).Y{w) = 

Vu^w : ^ Y(v!) - degree(u).F(it) = 

Y(v) = 1 

A non-negative set of values satisfying the above set is feasible for the dual LP. We find these by considering 
the resistive circuit S corresponding to the graph S (with each edge having unit resistance). The sink node 
is assigned ground potential and just enough current is injected at node w so that it attains unit potential. 
In terms of the standard potential functions described earlier, the potential that develops on any node u 
is s tt w (u). The potential value at node v, ix w (v), can be used to scale the input current at w thereby 
scaling all the potentials as well, such that potential at node v becomes unit. It follows that the values 
Y{u) = tt w (u)/tt w (v) and Y' equaling the value of the current injected form a feasible solution of the dual 
LP. The objective value at this point is, 



Y' H -r^ J~] h(u).TT w (u) 

TT W {V) ^ 

Since the current, Y', is a positive quantity, the above derivation implies the Lemma. □ 
We will now show that the lower bound is at most a constant factor smaller then the upper bound. We 
start with defining an auxiliary set of variables h v such that, 

h v = 1 = h v > 1 

From the positivity of 7r(.) function and the domination relation h v < h v , the following inequality follows. 



f= — \-^y^ u h u .-K W (u) < — — y^h a .7r w (u) (ii) 

v 7 u v ' u 

We will now bound the value F from below. 

Lemma 6.3. (Dhar Babai and Gorodezky J3j) In any stable recurrent configuration, for every edge, 

both the incident vertices cannot have zero particles. 

This follows from the fact that in any stable recurrent configuration, the last toppling of one of them 
would have taken place after the other and so the second node necessarily has at least one particle. For every 
edge u, v, at least one of h u and h v is > 1. 

Corollary 6.1. The set of nodes I = {v : h v = 0} form an independent set. 

Assume further that the graph satisfies (A), i.e. the maximum degree is A. Using this we will obtain the 
required bounds. Consider any vertex v € Vh along with its neighborhood (see figure ((6])). The function 7r(.) 
is harmonic over this neighborhood, so we have 

degree(v)7r(w) = ir(u) 
from which we get the sum of 7r(.) over any neighborhood in Vh as, 



ir(N(v)) = tt(v) + tt(u) = (degree(w) + l)ir(v) 

The contribution of any local region in the restricted potential trace T depends on just two possibilities 
regarding h v . We deal with them separately. 
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Case h v = 0: v £ I, consequently none of its neighbors are in /. In this case the contribution to T is 

the sum Y,u~v 7r («) = degree(u)7r(w) = a ^^+i 7r ( iV ( tl ))- 

Case /i„ = 1: v ^ I, and some of its neighbors are in /. Assume the worst case scenario when all the 
neighbors are in I, The contribution to T in such a situation is just ir(v) = c i cgrC e( 1 ,)-|-i 7r (^( t '))- 





Center node selected Center node's neighbors selected 

Figure 6: Neighborhood around vertex v; two cases, v selected and not selected in / 

It follows that any neighborhood N(v) contributes at least dcoTC e(. t ,) +1 tt(N(v)) to T, regardless of the 
specific values of h v . Since the degree is bounded by A, this translates to a minimum contribution of 
■^jtt(N(v)). Therefore, the value of T in equation (fTTj) . for any independent set /, is at least 



1 



(A + 1)it w (v) 

The upper bound yielded by Lemma |6. II is bounded from above by, 



(12) 



(A-l) 
ir w (v) 
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Using equations (fTU|) and (fT2")l . we obtain the following two sided bounds for degree bounded graphs, 



(A + iKAO ^ nwiu) ( ltw ^ Mu) (14) 

Theorem 6.1. For any sandpile with bounded vertex degrees, the minimum number of particles that need 
to be added at any vertex v to observe a toppling at any vertex w is equal, up to constant factors, to the 
following expression, 

j^jE'-M (15) 

The computation of transience class tcl(S') requires evaluating the above expression for all possible combi- 
nations of v, w € V(S). Computing the function 7r(.) can be done in very efficiently following the recent path 
breaking work by [28] , [29], [39J on solving symmetric, diagonally, dominant linear systems. Consequently, 
finding the pair with worst estimates is also easy to do. 

Corollary 6.2. There exists a polynomial time algorithm which computes the transience class of a degree 
bound sandpile up to constant factors. 

This answers one of the main questions raised in [10] . In this paper, the authors also derived several 
sandpile properties in analogy with the electric network theory. However, because of one sided bounds, most 
of these were also one sided. The sandwich bounds we have derived in this section will imply tighter versions 
for several of those potential theoretic properties for sandpile. 
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6.1 Proof of theorem 14.11 

We will now show that the local superposition property (LS) follows from mean value property (MV) and 
degree boundedness (A) of graph. 

Let H be the maximum number of particles which can be added at site v without causing a toppling at 
site w. Following the Theorem 16.11 and Lemma \6. II for degree bounded sandpile, the value of H is equal to 
(up to constant factors) 

Further, let h be maximum number of particles which when placed at each site of the ball B(v, r) do not 
topple w. In strict analogy with the discussion in Lemma 16.11 the following LP relaxation forms an upper 
bound on the value of h. 



max h 

Vm G B(v,r) : < z (u') — degree(it).z(u) + h < degree(u) — 1 

u' ™u 

Vu ^ B(v,r) : < z ( u ') — degree (w).z( it) < degree(it) — 1 

z(w) < 0,z > 0,h > 
The following minimization program is the dual of the above. 

min (degree (u) — l).Y(u) 

u 

Y(u ) + Y' - degree(u;).r(u;) > 

u' 

Mu^w: Y] Y{u') - degree(u).r(u) > 

E ^ 1 

u£B(v,r) 

Y > 0, Y' > 

We follow essentially the same argument as in Lemma 16.11 except that in the last step, we scale by 
J2 U £B(v r) 7r ( u ) instead of just ir(v). This leads to the following bound. 

h < £(d(u) - l).n w (u) (17) 

Lt,eB(»,r) KwW u 

We are given that the sandpile graph satisfies the mean value (MV) property (definition (jTU)0 . Which 
implies that, 

J2 n w (u)>C.n w (v).Vol(B(v,r)) (18) 

for some constant C . Using equations [17] and [HI we get the following bounds on h (up to constant 
factors) . 



ir w (v).Vol(B(v,r)) ^ 



(19) 
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Again, using the equation [TH] which tightly approximates the value of H and equation [T5] which bounds 
the value of h, we can bound h in terms of H. The following bound holds (up to constant factors). 



h * Vo K B( v ,r)) (20) 

We have thus demonstrated that the local superposition property (LS) follows from (MV) and (A). This 
completes the proof of theorem 14. II 



7 Special case of grid sandpile 

Figure (fTJ) shows an example of an m x n grid based sandpile. We consider the family of sandpiles consisting 
of the symmetric case (m = n). 

Consider an n x n grid graph. Attach an extra sink node s to all the nodes lying on the boundary of 
this grid such that there is a double edge with every corner node and single edge for others. We denote the 
sandpile and graph by GRID„. Figure {7J shows an example. 



(1,n) 



(n,n) 



(1,1) 




(n.1) 



Figure 7: An n x n grid sandpile : GRID n 



As in the example of grid sandpile, we assume that our sandpile is coming from the induced subgraph over 
a suitable set of vertices in the infinite lattice. The reason behind this choice is that we will be talking about 
the family of all grid sandpiles. While a particular finite grid can be a subgraph of many infinite graphs, 
the infinite grid is the canonical choice of a common ancestor of all finite grid graphs. We will show in this 
section that the infinite grid, Z 2 , along with the finite subgraph GRID„ satisfy all the properties required to 
imply polynomial bounds on the transience class of grids. 

Degree bound (A): The infinite grid Z 2 is regular with degree 4. Every induced subgraph obeys this degree 
bound. 

Polynomial volume growth property(V" a ): For any vertex p and radius r, volume of the ball B(p,r) 
grows as a quadratic in r. That is to say Z 2 satisfies V a with a = 2. This simple fact can be proved using 
elementary counting arguments. Figure ([7]) shows a typical ball around vertex p. 

Non-empty interior property (NI): Given any finite grid, GRID„, it suffices to show that there exists a 
path from center to every other node on which the distance to sink varies linearly. For simplicity, assume 
that n is odd and take any point v on the boundary. Denote the center node by c. Consider the canonical 
embedding of this grid in the plane. Join c to v by a line segment. Take the projection of this line segment 
towards x— axis and construct the path using the set of highest lattice points lying below the segment. This 
path has the requisite property and the growth parameter is just the slope of the line joining c and v. Figure 
([7]) shows one such path. The argument for the more general case of v lying inside the grid is identical. The 
case of even n can be handled using essentially the same ideas. 

This leaves us with showing the high local-conductance (hLC) and mean- value property (MV). The 
mean-value property is used in the Theorem I4.lt where it along with the property A, implies the Local 
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Superposition (LS) property. This is used in the proof of the Lemma B~T1 to prove the Overlapping Potentials 
(OP) property. Babai and Gorodezky [3] prove both (hLC) and (LS) property on grids, using elementary 
combinatorial arguments. We reproduce the discussion below without any essential changes. 

Definition 16. D4 symmetry: A function / defined on an infinite two dimensional lattice, Z 2 , / : Z 2 —> N 
is said to have D4 symmetry with respect to vertex v if is symmetric with respect to all the axis of symmetry 
passing through v. 

In case of functions which are defined on a finite n x n grid, the above definition is applied by assigning 
zero to all the points on which the function is not defined. This creates a function with a finite support. The 
point v, in these cases, is bound to be the geometric center of this finite grid. 

Definition 17. Axis Monotonicity: Let / be a function defined on an infinite two dimensional lattice, Z 2 , 
/ : Z 2 — > N. We say / is axis monotone about some vertex v if f(p) < f(q) for any pair of lattice points p 
and q such that the segment q — p is aligned perpendicular to some axis of symmetry passing through v and 
q is closer to this axis then p. 



See figure (|8} for a illustration of directions of the axis-monotonicity and D4 symmetry. 

Lemma 7.1. J^j/ If the starting (possibly unstable) configuration has D4 symmetry and is axis-monotone, 
then so is the resulting stable configuration. 

In the present context, the function of interest will be the number of particles present at a particular node. 
We observe that any function which is D4 symmetric, axis-monotone, and has finite support, is sandwiched 
inside an l\ ball (diamond) and l^ ball (square) of same radius. The number of sites contained in both these 
balls is quadratic in the radius. To be precise, denoting an l x ball of radius n around v by i?oo(i;,n), the 
volume of this ball is the maximum of particles it can hold in any stable configuration. This is given by, 



Similarly, the maximum weight of any stable configuration on B^v^n) is bounded by the following. 



The balls, therefore can hold at most 0(n 2 ) particles. When one starts adding particles at node v, clearly 
adding | B\(v,n) |= 0(n 2 ) particles is enough to flood an l\ ball of radius n around v. The property (hLC) 
follows immediately. 

The property (LS) requires more work. Here we start with placing some h particles at each node in 
B\(v,r). This configuration is axis- monotone and ^-symmetric. Consequently, the stable configuration 
resulting from this is also axis-monotone and £>4-symmetric. Let h be smallest such number such that the 
resulting configuration contains the ball Bi(v,R). The total number of particles this configuration can have 
is 0(i? 2 ). The particles we started of with are h.Q(r 2 ). Hence, h is <d((R/r) 2 ). The property (LS) follows. 




Figure 8: Axis Monotonicity of a function 



I BooM \=3(2n + l) 2 



I Bi(«,n) |= 6?i 2 + 6n + 3 
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Figure 9: A example of a function support which is ^-symmetric and axis monotone (about the center of 
the grid). The function support is shown in broken lines and bold lines depict the and l\ balls 

8 Parallels with classical complex analysis 

In the introductory section on potential theory on graphs, we briefly presented an example of multi-pole 
based construction of harmonic functions. Here we continue that discussion in the setting of real and complex 
analysis on continuous spaces. In analogy with the graphs induced on the vertex set Vh and the boundary 
sets connecting the poles to interior, we talk about open sets in R™ and the boundary of these sets. For the 
sake of illustration, consider an open set Q, C R 2 and its boundary set S£l. For definiteness, let f2 = B(p, r) 
that is a ball of radius r around the point p (Figure (fTOjl 'l. Consider a function h : 50. — > R + , that is an 
assignment of non-negative values to the boundary of this ball. Consider a real valued stochastic process 
{(Bi(t), -E?2(i)) : t > 0} with B\ and B2 independent and satisfying the following properties. 




Figure 10: A region fl inside R 2 : Brownian motion starting at the center p eventually hits the boundary <5f2 

• (£?i(0), ^2(0)) = p, the process starts at p. 

• For any pair of non-overlapping time intervals < t\ < ti < t% < £4, the increments, in both coordinates, 
Bifa) — Bi(t\) and B^{ti) — £^(£3) are independent. 

• For alii > and h > 0, Bi(t + h) — Bi(t) is normally distributed with mean zero and variance h. 

• The sample paths of this process are almost surely continuous. 

Above are the defining properties of a Brownian motion in two dimensions. We refer the reader to the 
recent comprehensive monograph by Morters and Peres |33) for a thorough discussion of the properties of 
Brownian motions and its analogies with random walks. Just like a random walk starting at some point hits 
any node on the boundary with some probability, a Brownian motion starting at some point hits any subset 
of the boundary with some probability. Using these probabilities, we can infer the distribution and, hence, 
the density of the measure induced by the Brownian motion on the boundary set. Note that the density is a 
function of the starting point of the Brownian motion. We now define the harmonic continuation of h over 
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the region f2. For every point q inside the ball, the value h(p) is defined as the expected value induced by 
this density on the boundary. 

A function h is called harmonic over Q iff it satisfies the Laplace equation (|2ip over all points inside fi. 




(21) 



One can also use the equivalent formulation of harmonic function in terms of its mean value property. 
Mean Value property of harmonic functions: If B(p,r) is the closed ball of radius r, centered at p, 
and completely inside £1 on which the harmonic function h is defined, then the value at the center of the ball 
is given by the average value at the surface of the ball. 



Here IS"" -1 ! is the total volume of the surface of the unit ball in n-dimensions and da is the n — 1 
dimensional surface measure. Conversely, all integrable functions satisfying the mean value property are 
harmonic. 

The equivalence of these conditions is discussed at length in [33]. In the context of graphs we refer the 
reader to the survey by Hoory, Linial and Wigderson [53] . They discuss and derive the concept of a Laplace 
operator in the context of graphs by first considering appropriate versions of divergence and gradient. Taking 
the divergence of the gradient, in the usual manner, yields the Laplacian. The linear operator corresponding 
to this is, not surprisingly, known as the combinatorial Laplacian of the graph. Once the formulations are in 
place, establishing the equivalence of mean value property and the Laplace equation is trivial. 

In our discussion till now, we have overlooked a very serious issue relating to boundary values as well as 
the boundary structure. The formulation in terms of Laplacian is local while that in terms of mean value 
is dependent on boundary values. The mean value formulation depends on the integral computation, which 
places geometric constraints on the structure of boundary sets. The property which the boundary is required 
to follow is called Poincare cone condition, (see [33], Definition 3.10). Once we ensure that the boundary 
set is well behaved, we also need a guarantee that for given boundary values, a harmonic continuation will 
exist. Indeed, this is the classical Dirichlet boundary value problem. Fortunately, the only condition required 
is that the values are continuous. The uniqueness of solution is easy to show. In fact the set of harmonic 
functions over a domain form a vector space. For a full discussion we refer the reader to the classical text by 
Rudin |35| . The most important examples of harmonic functions come from the real and imaginary parts of 
holomorphic (i.e. complex differentiable) functions. In this case, harmonicity of these functions is a simple 
consequence of the Cauchy-Riemann equations which every holomorphic function must satisfy. 

Following above discussion, it is reasonable to expect that a harmonic function should satisfy some non- 
trivial global properties. One such property is elucidated by the celebrated Harnack's inequality. We present 
below the two dimensional version for balls. 

Harnack's inequality for balls: If h is positive harmonic function defined over f2 then the following in- 
equality holds for any (open) ball B(p,r) C fi and q £ B(p,r). 



The generalization to higher dimensions is straightforward. From this, the inequality for all compact 
connected regions K(C. f2) follows using a very simple chaining argument which involves considering a set of 
open balls in K which cover the path from the reference vertex p to the vertex q. Using the bounds from 
equation ([23| . the growth rate of harmonic functions in each of these balls is bounded by the respective 
constants Hi. Multiplying these gives us an obvious bound over the considered path. By taking any finite 
covering of K (compactness of K guarantees the existence of such a covering) , one can find an absolute bound 
on the value of the Harnack's constant by simply multiplying the constants corresponding to each of the balls 
(see figure This implies the Harnack's inequality over any compact region. 

Harnack's inequality over general regions: Suppose that f2 is connected and if is a compact subset. Then 
for any harmonic function, h, defined over f2, the following inequality holds 




(22) 




(23) 
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K 




h(q) < H K h{p) (24) 

Here Hx is the Harnack's constant of the region K and is independent of the function h. Note that the 
version for balls (|23p uses the exact value for this constant. This value comes from an alternate formulation of 
the integral equation (f2"2"| based on the Poisson kernel. For further discussion on the Poisson kernel, we refer 
the reader to Rudin |35j . Harnack's inequality plays fundamental role in the traditional complex analysis 
and PDE theory. Recently it found use in Perelman's proof of the Poincare conjecture. 

The Harnack's constant, in a sense, summarizes the connectivity of the corresponding region. To make 
sense of this statement, we consider the analogous situation in graphs. Consider the infinite grid graph, 
denoted by Z 2 . Consider any finite subset of vertices Vh and the subgraph induced on them, in the sense 
discussed in previous sections. See figure (fT2"j) . The black nodes belong to Vh, heavy edges are incident on 
sink and for clarity the selected sub-grid is demarcated by the light thick line. The Harnack's constant for 
this region should measure the maximum possible growth in magnitude of a function which is harmonic over 
this domain. This is equivalent to finding the maximum value of 7r p (g) - 1 , with the sink node as the zero 
potential, and both p and q on the vertex boundary. Consider one such pair p and q (as shown in figure 
(|12p ). The value ^ P (q) measures the probability of a walk, starting at p, finding q before it dissipates out 
of the region. A high value (bounded from below by some constant, say) would imply that the nodes under 
consideration enjoy overall good connectivity and even a memory less strategy can find some path between 
them. A high value of 7r p (g) , on the the other hand, would imply that the pair is not so well connected in 
the potential theoretic sense and it is highly likely for a walk starting at p to exit the system before hitting q 
at all. Bounding the value of 7r p (g) _1 for all possibilities of p and q implies that even the pair with the worst 
case connectivity are not very bad. These properties are analogous to the mixing times of graphs without 
boundaries. 

Remark: The Figure (jTTJ) shows the path from p to q covered by a set of compact balls. One can see that even 
if these balls have significant overlapping, the discussion in preceding paragraphs is oblivious to the utility. 
One might as well have a set of balls, any pair of which intersects at most at one point. Evidently, the point 
to point connectivity, Brownian motion based diffusion type, is very different for the two cases and when 
the covering balls have significant overlapping, the Harnack's constant must be (much) smaller. It therefore 
seems desirable to formulate properties which would allow efficient gluing of Harmonic functions. For the 
continuous case, no such properties are known. Recall the triangle inequality for potentials (Lemma 13. 2j) on 
graphs. The analogous version for continuous spaces is almost as old as the classical potential theory itself. 
This is enough to perform the stringing of Harnack's constants mentioned above. In the context of graphs, 
when we try to bound the transience classes, relying on this property alone clearly doesn't suffice. One can 
show that, in case of an n x n grid, about O(logn) balls will be required to cover the path between opposite 
corners. Multiplying their transience classes yields super-polynomial bounds. The impulse superposition 
properties we derive, in the context of sandpile, are essential to show efficient bounds and form a completely 
new component in any effort towards demonstrating gluing properties of potential functions. Proving these 
properties in the context of graphs is more general then the traditional setting of continuous spaces. This 
is so because there exists the canonical lattice graph for every Euclidean space and the analogies between 
the properties of these two objects are well known. However, there exist many graphs, for example various 
fractal graphs, that do not arise as discretization of any continuous space. 
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Figure 12: A finite subset of the grid graph 

In the preceding sections, we demonstrate polynomial bounds on the transience class when the graph sat- 
isfies certain assumptions. We now discuss these assumptions, and try justifying these by showing analogous 
constraints in the continuous framework of harmonic analysis. 

- Volume growth parameter, a, ensures uniformity in measure distribution throughout the space. Its 
role is analogous to the dimension of spaces in classical analysis. Whether we consider the counting 
measure associated with graphs and or the Lebasgue measure associated with continuous spaces, volume 
growths in both the cases are polynomials with exponents a, and is uniform throughout the space under 
consideration. In the derivation of our bound on the transience class and consequently on the Harnack's 
constant of the graph, this analogy is evident. 

- Mean Value Property (MV) of graphs is similar to the mean value property of continuous spaces. This 
assumption underlines a certain homogeneity property of space and uniformity in measure distribution. 
The main implication of such homogeneity is that a random walk which starts at center of a ball, does 
not significantly favor any point (or set of points) on the surface over others. 

- Non-Empty Interior Property (NI) is an analogue of the assumption of non-degeneracy of domains in 
analysis. Which means that for every point inside the region, there is a full-dimensional ball containing 
it which lies completely inside the region. The concluding section also discusses degeneracy in graphs. 

- High local Conductance (hLC) is analogous to ensuring that the measure distribution in space is 
(kind of) isotropic. This means that when diffusion starts from any point, for any potential value, 
the distances of nearest and farthest equipotential points from the center are not much different. In 
continuous spaces, this is a simple consequence of symmetry. In fact, one can show much stronger 
property of conformal invariance, in the context of Brownian motions. 

- The Polynomial Transience Property (pT) is similar to ensuring that the Harnack's constant is not 
too large. Following from our previous discussion, this implies that the diffusion inside the region is 
reasonably good. To be precise, polynomial bounds on transience class mean that the corresponding 
sandpile achieves recurrence in reasonable time, and therefore is a good model of diffusion. Conversely, 
a sandpile with exponential transience class should correspond to a system which has not been modeled 
well. We will continue this discussion in the following sections. 
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9 Future work 



In the present work, we generalize the results of Babai and Gorodezky [3J beyond the family of grid graphs. 
The analysis raises lot of interesting questions, answers to which will aid our understanding of discrete 
diffusions. 

The first question is, how far are we from a complete characterization? The present result shows that 
under some very general and reasonable assumptions, the polynomial bounds can be established. There 
are some obvious extensions possible from our theorems. For example if one considers graphs made from 
attaching two graphs each individually satisfying the required pre-conditions but having different volume 
growth parameters. One can bound the Harnack's constant for both components separately and bound the 
one for the whole graph using the simple triangle inequality. But apart from these trivial modifications, we 
do not yet know of any class of sandpile which are not covered by our result. Supported by the analogy with 
the continuous case, we believe the assumptions we make are essential and can not be substantially relaxed. 

The second question supplements the first one. Consider any Euclidean space and triangulate it using 
simplices or cubic complexes, basically any regular tessellation will do. We are interested in the graph of 
this triangulation complex. Indeed, the grid sandpile, which started of the study of transience class bounds, 
are obtained in exactly this manner from the two dimensional Euclidean space. In this graph, consider a 
subgraph corresponding to any open set. To be precise, the mesh size is also a parameter. So for every set 
in space, there exists a family of subgraphs parameterized by the fine-ness of triangulation. For example, 
in the case of grid sandpile, the underlying set is a square and the parameter n in the n x n grid sandpile 
encapsulates the precision with which the discrete version can model the continuous version. In this setting, 
one can see that polynomial bounds in sandpile settings correspond to non-trivial bounds on the Harnack's 
constant in the continuous space. Any degenerate set in K™ has unbounded Harnack's constant. It can be 
shown that the sandpile corresponding to these sets also have transience classes which cannot be polynomially 
bounded. Babai and Gorodezky [3] first mention this for line sandpile. The case of general degenerate sets 
is similar. The only crucial component is the fact that from every node, there is an edge to sink node. One 
can see an interesting confluence of the notion of bad models here. In the setting of topological spaces, full 
dimensional sets are considered to be available or good models (physically speaking) while the degenerate ones 
are considered beneath consideration as they are not good models (in whichever context). On similar lines, 
Babai and Toumpakari [5], and Babai and Gorodezky [3J make case for the notion of polynomially bound 
transience classes as the essential qualification for a sandpile to be a good model. Dhar [TB] mentions several 
physical phenomena which sandpiles are intended to model. In all these settings, the sandpile graphs arise 
as tessellation graphs of the underlying continuous space. It is only natural to expect that if the underlying 
set itself is a bad model, the sandpile one top of it is going to a bad model. 

The line sandpile is a simplest case of this observation. This sandpile is based on a line graph which has 
two connections to sink from every node, except the corner ones which have three connections each. One can 
think of these sandpiles as n x 1 grids based sandpiles (or degenerate grids, to make the connection clearer). 
Babai, Gorodezky, and Shapiro [3] show exponential lower bounds on the transience classes of these sandpiles 
using very simple combinatorial arguments. Independently, Choure and Vishwanathan |10) also demonstrate 
these same bounds, using potential theoretic arguments. With some work, these can be extended to other 
degenerate cases. This categorically classifies the line sandpile as an infeasible model for diffusion. Now 
consider the set in R 2 whose tessellation gives rise to such graph. The underlying set is a line segment 
in a plane, clearly a degenerate set. In this sense, line sandpile is a misnomer. The correct sandpile to 
associate with this name should be the one corresponding to a line segment in M. This one has just two single 
connections to sink node, one at each end, and not surprisingly, has polynomially bound transience class. In 
fact any sandpile family of k x n grids, where k is a constant, cannot correspond to an open set and should 
have exponential transience classes. Empirical results support this view. This leads us to a meta-conjecture 
of sorts that impinges on the notion of a dimension theory for graphs itself. 
Question: Does every bad sandpile correspond to a degenerate set? 

The constant factor approximation which we discuss in this paper, demonstrates that studying the tran- 
sience class is similar to bounding the harmonic functions over the underlying graphs. The growth rates of 
harmonic functions in continuous spaces provide important information about the dimension of the space, as 
already discussed in the previous section. One can generalize this notion to graphs. The main problem being 
that unlike Euclidean spaces, graphs can be highly non-uniform. So we restrict our attention to graphs which 
are uniform in terms of measure distribution. Imposing same volume growth laws at every vertex ensures 
this. There are graphs which do not fall in this category but locally, in well defined regions, satisfy the 
uniformity we seek. In this case one can analyze the regions separately and then try to deduce the properties 
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of the whole graph from these components. Finally we have highly heterogenous graphs. Here, the geometric 
notion of dimension completely breaks down. Since these graphs arise in completely combinatorial settings, 
properties like dimension are much harder to attribute and justify. Coming back to the measure uniform 
graphs, the question is how much geometric information can be conveyed through the transience class? 

The dimension of a graph and deciding degeneracy of subgraphs are obvious applications. The theory 
of Ricci flows on Reimannian manifolds has recently found immense favor in mathematical community. To 
develop the analogous theory on graphs, the Harnack's constant is potentially going to play significant role. 
The second related question in the study of harmonic functions is the gluing property we discuss in last 
section. Currently there is no known result in this direction. While we obtain a decent bound, the proof goes 
though sandpiles and a lot of apparently unrelated combinatorics. Is there a way to derive a gluing theorem 
for harmonic functions, even in just two dimensional space without resorting to sandpiles? 

To conclude, we bring together ideas from three different areas of mathematics. Our bounds demonstrate 
that the research in each of these areas has ramifications in the other two. In this sense, the importance of 
resolving any question in any of these areas is greatly magnified. We hope that this paper in conjunction 
with our earlier work motivates further investigation on discrete diffusions, in analogy with random walks 
and complex analysis. 
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